//
//  Simulation.cpp
//  phy_four
//
//  Created by Snir Gazit on 3/12/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#include "Simulation.h"

typedef boost::chrono::system_clock::time_point TimePoint;
typedef boost::chrono::duration<double, boost::ratio<3600>> Hour;
typedef boost::chrono::duration<double, boost::micro> MicroSec;

double calcWR(double mu, double U, int m) {
	double c1, c2, c3, md = double(m), arg;
	arg = 0.25 * pow(mu, 2) / U;
	c1 = 0.25 * pow(U, -(5. / 4.) - double(md) / 2);
	c2 = mu * gsl_sf_gamma((5. / 4.) + md / 2)
			* gsl_sf_hyperg_1F1((5. / 4.) + md / 2, 3. / 2., arg);
	c3 = sqrt(U) * gsl_sf_gamma((3. / 4.) + md / 2)
			* gsl_sf_hyperg_1F1((3. / 4.) + md / 2, 1. / 2., arg);
	return c1 * (c2 + c3);
}

void calcWeight(double muIn, double UIn, double g, int maxBond,
		ostream &outFile, BoltzWeightType &Rweight) {
	double mu = muIn * g, U = UIn * g;
	double n, c, w;
	c = calcWR(mu, U, 0);
	for (int i = 1; i < MAX_BOND; i++) {
		n = calcWR(mu, U, i);
		w = (n / c);
		c = n;
		Rweight[i - 1] = w;
	}
}
// Worm Alg

void OnModel::init(QmcParameters &qmcParameters) {
	//Init Params
	qmcParameters_ = qmcParameters;
	string runNum = qmcParameters_.runNum.getVal();
	// OutFile
	String fileName = qmcParameters_.outPath.getVal() + "/out" + runNum
			+ ".txt";
	outFile_.open(fileName.c_str());
	outFile_ << "Starting N=3 Simulation, Run Number - " + runNum;
	BoltzWeightType RWeight;
	// Calc Weights
	qmcParameters_.g.setVal(1. / (qmcParameters_.g.getVal()));
	if (qmcParameters_.simType.getVal() == "HEISENBERG")
		RWeight.fill(1.);
	else
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(),
				qmcParameters_.g.getVal(), qmcParameters_.maximalBonds.getVal(),
				outFile_, RWeight);
	double switchProb = 0.5;
	double advanceProb = (1 - switchProb) / 2.;
	boost::array<double, 3> prob = { switchProb, advanceProb, advanceProb };
	openClosedSwitch_.init(&eng_, DiscreteVariableParams(prob));
	// Init Lattice
	boost::array<int, DIM> sizes;
	sizes.fill(qmcParameters_.L.getVal());
	if (qmcParameters_.simType.getVal() == "HIGGS")
		sizes[DIM - 1] = qmcParameters_.T.getVal();
	lat_.init(sizes);
	lat_.setSectorType(GXY);
	lat_.getIra()->XYTotalBonds_++;
	// Init Rng
	eng_.seed(qmcParameters_.seed.getVal());
	// Init Moves
	//	createMoveXY_.init(&eng_, &lat_, qmcParameters_, XYWeight);
	//	createMoveIsing_.init(&eng_, &lat_, qmcParameters_, IsingWeight);
	//  removeMoveXY_.init(&eng_, &lat_, qmcParameters_, XYWeight);
	//removeMoveIsing_.init(&eng_, &lat_, qmcParameters_, IsingWeight);
	switchMoveIsing_.init(&eng_, &lat_, qmcParameters_, RWeight);
	switchMoveXY_.init(&eng_, &lat_, qmcParameters_, RWeight);
	shiftMoveIsing_.init(&eng_, &lat_, qmcParameters_, RWeight);
	shiftMoveXY_.init(&eng_, &lat_, qmcParameters_, RWeight);
	jumpMoveXY_.init(&eng_, &lat_, qmcParameters_, RWeight);
	jumpMoveIsing_.init(&eng_, &lat_, qmcParameters_, RWeight);
	chooseLoopType_.init(&eng_, BinomialVariableParams(0.5));
	// Measurements
	measures_.init(&lat_, qmcParameters_, &RWeight, &outFile_);
	// Read From File
	if (qmcParameters_.initFromFile.getVal())
		this->read();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	//measures_.updateMinNumOfMeasure();
}
OnModel::~OnModel() {
	outFile_.close();
}

void OnModel::thermalize() {
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	acc CProb, XYProb;
	for (int i = 0; i < numOfThermal; i++) {
		nextMove();
		if (lat_.isClosed())
			CProb(1.);
		else
			CProb(0.);
		if (lat_.getSectorType() == GXY)
			XYProb(1.);
		else
			XYProb(0.);
	}
	bool conv;
	double cCProb = ba::mean(CProb), nCProb = 0.0;
	do {
		CProb = acc();
		XYProb = acc();
		for (int i = 0; i < numOfThermal; i++) {
			nextMove();
			if (lat_.isClosed())
				CProb(1.);
			else
				CProb(0.);
			if (lat_.getSectorType() == GXY)
				XYProb(1.);
			else
				XYProb(0.);
		}
		double nCProb = ba::mean(CProb);
		outFile_ << "Closed Prob " << nCProb << endl;
		outFile_ << "XY Prob " << ba::mean(XYProb) << endl;
		double ratio = fabs(1. - (cCProb / nCProb));
		conv = ratio < 0.01;
		cCProb = nCProb;
	} while (!(conv));
	//CreateLoopMove_.setAvWormLength(2 * ba::mean(wormLength));
	outFile_ << "Final Closed Prob - " << cCProb << endl;
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) * cCProb);
	//	closedRate = 1;
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}

bool OnModel::run() {
	bool endTime, dumpFlag;
	namespace ba = boost::accumulators;
	ba::accumulator_set<double, ba::stats<ba::tag::mean> > wormLength;
	while (measures_.checkNotDone()) {
		nextMove();
		measures_.DoMeasurement(endTime, dumpFlag);
		if (dumpFlag) {
			this->dump(false);
		}
		if (endTime) {
			//this->dump(true);
			break;
		}
	}
	return !endTime;
}
void OnModel::nextMove() {
	SectorType cSectorType = lat_.getSectorType();
	if (cSectorType == GISING) {
		// Check in Closed
		if (lat_.isClosed()) {
			// Jump Or Shift
			int openClosedSwitch = openClosedSwitch_.draw();
			if (openClosedSwitch == 0)
				switchMoveIsing_.makeMove();
			else if (openClosedSwitch == 1)
				jumpMoveIsing_.makeMove();
			else
				shiftMoveIsing_.makeMove();
			// Shift
		} else {
			shiftMoveIsing_.makeMove();
		}
	} else {
		// Check in Closed
		if (lat_.isClosed()) {
			// Jump Or Shift
			int openClosedSwitch = openClosedSwitch_.draw();
			if (openClosedSwitch == 0)
				switchMoveXY_.makeMove();
			else if (openClosedSwitch == 1)
				jumpMoveXY_.makeMove();
			else
				shiftMoveXY_.makeMove();
			// Shift
		} else {
			shiftMoveXY_.makeMove();
		}

	}
}
void OnModel::dump(bool final) {
	String outPath = qmcParameters_.outPath.getVal();
	measures_.dump(outPath, final);
	lat_.dump(outPath, final, &outFile_);
	eng_.dump(outPath);
}
void OnModel::read() {
	String outPath = qmcParameters_.outPath.getVal();
	measures_.read(outPath);
	lat_.read(outPath);
	eng_.read(outPath);
}

// WormLoop Algorithm

// Worm Alg
void OnModelLoop::init(QmcParameters &qmcParameters) {
	//Init Params
	qmcParameters_ = qmcParameters;
	string runNum = qmcParameters.runNum.getVal();
	// OutFile
	String fileName = qmcParameters_.outPath.getVal() + "/out" + runNum
			+ ".txt";
	outFile_.open(fileName.c_str());
	outFile_ << "Stating N=3 Simulation - Run Number - " + runNum << endl;
	if (qmcParameters_.simType.getVal() == "XY") {
		//boltzRatio_ = BoltzWeightType(qmcParameters_.maximalBonds.getVal(), 1.);
		boltzWeight_.fill(1.);
		qmcParameters_.g.setVal(1 / (2. * qmcParameters_.g.getVal()));
	} else {
		qmcParameters_.g.setVal(1. / (qmcParameters_.g.getVal()));
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(),
				qmcParameters_.g.getVal(), qmcParameters_.maximalBonds.getVal(),
				outFile_, boltzWeight_);
	}
	// Init Lattice
	lat_.init(qmcParameters_.Llist.getVal());
	// Init Rng
	eng_.seed(qmcParameters_.seed.getVal());
	// Init Moves
	CreateLoopMoveXY_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	CreateLoopMoveIsing_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	chooseLoopType_.init(&eng_, BinomialVariableParams(1, 1. / 2.));
// Measurements
	measures_.init(&lat_, qmcParameters_, &boltzWeight_, &outFile_);
// Read From File
	if (qmcParameters_.initFromFile.getVal())
		this->read();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	measures_.updateMinNumOfMeasure();
}
void OnModelLoop::setCoupeling(double g) {
	if (qmcParameters_.simType.getVal() == "HEISENBERG")
		boltzWeight_.fill(1.);
	else
		calcWeight(qmcParameters_.mu.getVal(), qmcParameters_.U.getVal(), g,
				qmcParameters_.maximalBonds.getVal(), outFile_, boltzWeight_);
	CreateLoopMoveIsing_.setCoupeling(boltzWeight_, g);
	CreateLoopMoveXY_.setCoupeling(boltzWeight_, g);
}
OnModelLoop::~OnModelLoop() {
	outFile_.close();
}

void OnModelLoop::thermalize() {
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	for (int i = 0; i < 10 * numOfThermal; i++)
		nextMove();
	double cLengthXY = ba::mean(wormLengthXY_), nLengthXY;
	double cLengthIsing = ba::mean(wormLengthIsing_), nLengthIsing;
	double cTotal = cLengthXY + cLengthIsing, nTotal, cmaXY, cmaIsing, nmaXY,
			nmaIsing;
	bool convXY, convIsing, convRatio;
	acc maXY, maIsing;
	maXY(cLengthXY);
	maIsing(cLengthIsing);
	cmaXY = ba::mean(maXY);
	cmaIsing = ba::mean(maIsing);
	do {
		wormLengthIsing_ = acc();
		for (int i = 0; i < numOfThermal; i++)
			nextMove();
		nLengthXY = ba::mean(wormLengthXY_);
		nLengthIsing = ba::mean(wormLengthIsing_);
		maXY(nLengthXY);
		maIsing(nLengthIsing);
		nTotal = nLengthXY + nLengthIsing;
		nmaXY = ba::mean(maXY);
		nmaIsing = ba::mean(maIsing);
		outFile_ << "Worm Length XY MA - " << nmaXY << " Worm Length Ising MA"
				<< nmaIsing << " Total MA - " << nmaXY + nmaIsing << endl;
		double ratioXY = fabs(1. - (nmaXY / cmaXY));
		double ratioIsing = fabs(1. - (nmaIsing / nmaIsing));
		double ratioBoth = fabs(1. - (nmaXY / nmaIsing));
		convXY = ratioXY < 0.05;
		convIsing = ratioIsing < 0.05;
		convRatio = ratioBoth < 0.1;
		cmaIsing = nmaIsing;
		cmaXY = nmaXY;
	} while (!(convXY && convIsing && convRatio));
	outFile_ << "Final Worm Length XY - " << ba::mean(wormLengthXY_)
			<< " Final Worm Length Ising " << ba::mean(wormLengthIsing_)
			<< " Final Worm Length Total "
			<< ba::mean(wormLengthXY_) + ba::mean(wormLengthIsing_) << endl;
	double wormLength = nTotal;
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) / wormLength);
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}
void OnModelLoop::simulatedAnnealing() {
	outFile_ << " Startting Simulated Anneling " << endl;
	double gHighT = 1. / qmcParameters_.gHighT.getVal();
	int numOfSMA = qmcParameters_.numOfGSA.getVal();
	double dg = (qmcParameters_.g.getVal() - gHighT) / double(numOfSMA - 1);
	double nLengthIsing, nLengthXY, nTotal;
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	vector<double> gSA;
	for (int i = 0; i < qmcParameters_.numOfGSA.getVal(); i++)
		gSA.push_back(gHighT + dg * double(i));
	for (auto it_g = gSA.begin(); it_g != gSA.end();) {
		double cg = (*it_g);
		setCoupeling(cg);
		outFile_ << "Current g -" << cg << endl;
		bool conv;
		int numIter = 0;
		do {
			wormLengthXY_ = acc();
			wormLengthIsing_ = acc();
			for (int i = 0; i < numOfThermal; i++)
				nextMove();
			nLengthXY = ba::mean(wormLengthXY_);
			nLengthIsing = ba::mean(wormLengthIsing_);
			nTotal = nLengthXY + nLengthIsing;
			numOfThermal = int(
					100 * double(lat_.getNumOfSites()) / double(nTotal));
			outFile_ << "Worm Length XY MA - " << nLengthXY
					<< " Worm Length Ising MA - " << nLengthIsing
					<< " Total MA - " << nTotal << " numOfThermal "
					<< numOfThermal << " Iter " << numIter << endl;
			double ratio = (nLengthXY / nLengthIsing);
			conv = ((ratio < 1.5) && (ratio > 0.75));
			numIter++;
		} while ((!conv)); //&& (numIter < 50));
//		if (conv)
		it_g++;
		//	else {
		//	if (it_g != gSA.begin())
		//	it_g--;
		//	}
	}
	numOfThermal = int(5000 * double(lat_.getNumOfSites()) / double(nTotal));
	wormLengthXY_ = acc();
	wormLengthIsing_ = acc();
	for (int i = 0; i < numOfThermal; i++)
		nextMove();
	nLengthXY = ba::mean(wormLengthXY_);
	nLengthIsing = ba::mean(wormLengthIsing_);
	nTotal = nLengthXY + nLengthIsing;
	double wormLength = nTotal;
	int closedRate = int(
			double(qmcParameters_.rateFactor.getVal())
					* double(lat_.getNumOfSites()) / wormLength);
	measures_.updateClosedRate(closedRate);
	outFile_ << "Closed Rate - " << closedRate << endl;
}
bool OnModelLoop::run() {
	bool endTime, dumpFlag;
	namespace ba = boost::accumulators;
	wormLengthXY_ = acc();
	wormLengthIsing_ = acc();
	while (measures_.checkNotDone()) {
		nextMove();
		measures_.DoMeasurement(endTime, dumpFlag);
		if (dumpFlag) {
			this->dump(false);
		}
		if (endTime) {
			//this->dump(1);
			break;
		}
	}
	outFile_ << "Final Worm Length XY - " << ba::mean(wormLengthXY_)
			<< " Final Worm Length Ising " << ba::mean(wormLengthIsing_)
			<< " Final Worm Length Total "
			<< ba::mean(wormLengthXY_) + ba::mean(wormLengthIsing_) << endl;
	outFile_ << " Erase Prob - " << CreateLoopMoveXY_.getEraseProb() << endl;
	return !endTime;
}
void OnModelLoop::nextMove() {
	if (chooseLoopType_.draw() == 0)
		wormLengthXY_(CreateLoopMoveXY_.makeMove());
	else
		wormLengthIsing_(CreateLoopMoveIsing_.makeMove());

}
void OnModelLoop::dump(bool final) {
	String outPath = qmcParameters_.outPath.getVal();
	measures_.dump(outPath, final);
	lat_.dump(outPath, final, &outFile_);
	eng_.dump(outPath);
	String fileName = outPath + "/sim_params.xml";
// Write to XML
	pugi::xml_document doc;
	pugi::xml_node rootNode = doc.append_child("SIM_PARAMS");
// Parameters
	pugi::xml_node node = rootNode.append_child("CLOSED_RATE");
	node.append_child(pugi::node_pcdata).set_value(
			(boost::lexical_cast<string>(measures_.getClosedRate())).c_str());
	doc.save_file(fileName.c_str());
//	doc.save(std::cout);
//	cout << endl;
}

void OnModelLoop::read() {
	outFile_ << "Starting Reading From File" << endl;
	String outPath = qmcParameters_.outPath.getVal();
	measures_.read(outPath);
	outFile_ << "Measures Loaded" << endl;
	lat_.read(outPath);
	int totalBonds;
	try {
		bool testConf = lat_.testConfiguration(totalBonds);
		outFile_ << "Total Bonds Calculated - " << totalBonds
				<< " total Bonds previous sim - " << lat_.getTotalBonds()
				<< endl;
		if (totalBonds != lat_.getTotalBonds())
			throw "Total Bonds after Read don't match";
		if (!testConf)
			throw "Bad Load of Configuration";
	} catch (string &e) {
		outFile_ << e;
		exit(1);
	}
	outFile_ << "Lattice Loaded" << endl;
	eng_.read(outPath);
	String fileName = outPath + "/sim_params.xml";
// Read From XML
	pugi::xml_document doc;
	try {
		pugi::xml_parse_result result = doc.load_file(fileName.c_str());
		if (!result)
			throw result;
	} catch (pugi::xml_parse_result &result) {
		std::cout << "Error description: " << result.description() << "\n";
		exit(1);
	}
	pugi::xml_node rootNode = doc.child("SIM_PARAMS");
	measures_.updateClosedRate(readXmlData<int>(rootNode, "CLOSED_RATE"));

}
